arXivrl 506.08089V 1 [astro-ph.EP] 26Jun2015 


Realistic uncertainties on Hapke model parameters 
from photometric measurement 


Schmidt Frederic^Fernando Jennifer^ 

^ Univ. Paris-Sud, Laboratoire GEOPS, UMR 8 I 4 . 8 , Orsay, F-91405; 
(frederic.schmidt@u-psud.fr) ^ CNRS, Orsay, F-91405, France 


Abstract 

The single particle phase function describes the manner in which an average el¬ 
ement of a granular material diffuses the light in the angular space usually with 
two parameters: the asymmetry parameter b describing the width of the scat¬ 
tering lobe and the backscattering fraction c describing the main direction of the 
scattering lobe. Hapke proposed a convenient and widely used analytical model 
to describe the spectro-photometry of granular materials. Using a compilation 
of the published data, Hapke (2012, Icarus, 221, 1079-1083) recently studied 
the relationship of b and c for natural examples and proposed the hockey stick 
relation (excluding b > 0.5 and c > 0.5). For the moment, there is no theoretical 
explanation for this relationship. One goal of this article is to study a possible 
bias due to the retrieval method. 

We expand here an innovative Bayesian inversion method in order to study 
into detail the uncertainties of retrieved parameters. On Emission Phase Func¬ 
tion (EPF) data, we demonstrate that the uncertainties of the retrieved param¬ 
eters follow the same hockey stick relation, suggesting that this relation is due 
to the fact that b and c are coupled parameters in the Hapke model instead of 
a natural phenomena. Nevertheless, the data used in the Hapke (2012) com¬ 
pilation generally are full Bidirectional Reflectance Diffusion Function (BRDF) 
that are shown not to be subject to this artifact. 

Moreover, the Bayesian method is a good tool to test if the sampling geom¬ 
etry is sufficient to constrain the parameters (single scattering albedo, surface 
roughness, 5, c, opposition effect). We performed sensitivity tests by mimicking 
various surface scattering properties and various single image-like/disk resolved 
image, EPF-like and BRDF-like geometric sampling conditions. The second 
goal of this article is to estimate the favorable geometric conditions for an accu¬ 
rate estimation of photometric parameters in order to provide new constraints 
for future observation campaigns and instrumentations. 

Keywords: spectro-photometry, Hapke model, Bayesian inversion, disk 
resolved image, EPF, BRDF, uncertainties 
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1. Introduction 


Photometry is the study of the surface by the angular response of the re¬ 
flected light by a medium described by the Bi-directional Reflectance Distribu¬ 


tion Function (BRDF) (Hapke 1993) 


Hapke proposed a set of approximated analytical equations to estimate con- 


veniently the BRDF of a granular medium (e.g. Hapke 1981 Hapke and Wells 

1981 Hapke 19841 1986 2002 2008 

). This formulation has been controversial 

for two decades (e.g. Mishchenko 1994 ] 

Hapkel 1996 Shepard and Helfenstein 

|2007 Shkuratov et al. 2012 Hapke 

20R 

)) but due to its relative simplicity and 


fast computation, many authors have been using it to analyze laboratory data 


(e.g. 

Cord et al. 2003 Souchon et al. 2011 Beck et al. 

2012 

Pommerol et al. 

pMs 

Johnson et al. 2013), telescopic observation (e.g. 

Hapke et ahj |1998 

), in 

situ data (e.g. Johnson et al.l 1999| 2006b|a), remote sensing data (e.g. Jehl 

et al.[ |2008| lYokota et al.[ |2011| |Fernando et al.[ 2013| Vincendon[ |2013| Sato 


et ah 2014) 


The scope of this article is to discuss the properties of the Hapke model 
in term of data analysis, but not to discuss particular aspects of the realism 
of the photometric Hapke model. Some authors addressed the difficulties to 
fit the model to actual data, indicating that the problem is ill posed, due to 


parameter coupling (Mustard and Pieters 1989 Helfenstein and Veverka 1989 


Baratoux et ah 2006 Souchon et ah 2011). The most usual ways to fit are 


a minimization of the stepwise in a grid (Shepard and Helfenstein 


using the Levenberg-Marquardt method (|Sato et al. 2014 Johnson et al. 


a simplex algorithm (Gunderson et al. [2006 ), and a genetic algori thm ([Cord 


et al. 2003) and Particle Swarm Optimization (Beck et al. 2012 Pommerol 


2007) 


2013) 


et al.[ |2013[ ) . These strategies are relevant in the case of an unique solution (only 
one minimum in the y^) and close to gaussian shape around the solution. Due 
to the non-linearities of the Hapke equations, those two mathematical properties 
may not be fulfilled. 

We propose here a new kind of technique, based on the Bayesian formula¬ 


tion to estimate the model parameters in agreement with the data (Tarantola 
[and Valette , 1982). The general framework of Bayesian theory enabled us to: 


(i) take into account precisely the uncertainties of measured quantities, (ii) de¬ 
fine precisely the range of possibility of all model parameters, (hi) estimate the 
range of solution in a general case (that may not have a gaussian shape). From 
theoretical point of view, each information is described as a probability den¬ 
sity function (PDF): the measured quantities, the a priori parameters and the 
posterior parameters. Within this framework, the solution is expressed as a 
“final state of information” which always exists, solving the apparent ill-posed 
problem. 

We already used this strategy to analyze the spectro-photometric data from 
the Compact Reconnaissance Imaging Spectrometer for Mars (CRISM) instru¬ 
ment (Murchie et al. 2007), after the MARS-ReCO atmospherical correction 


(Ceamanos et al. 2013), by estimating the PDF of all Hapke parameters of 


the Martian surface (Fernando et al. 2013). We also applied it on full CRISM 
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images on the MER landing sites (Fernando et al. 2015) 


In this article, we propose to perform extensive sensitivity tests on synthet¬ 
ical dataset, using the Bayesian inversion, in order to study the behavior of 
Hapke model to fit the data. This study aims at: 


1. estimating precise uncertainties level on the model parameters on differ¬ 
ent typical observation types (one image-like observation/disk resolved 
image, CRISM-like Emission Phase Function (EPF), very favorable EPF, 
Bidirectional Reflectance Distribution Function (BRDF) measurements) 

2. determining the uncertainties dependance on BRDF sampling 

3. evaluating the possibility to have more than one “solution”, i.e. multiple 
minima 

4. chasing any effect which could bias the estimation of the parameters b and 
c leading to a fake hockey stick relationship 

This work should help to interpret previous analyses but also to design future 
instrumental and observational campaigns. 


2. Method 


2.1. Direct model: the Hapke^s photometric model 


Standard 1993 Hapke modeling (Hapke 1993) is widely used in the planetary 


science community due to the simplicity of its expression, fast computation, and 
the purported physical meaning of the model parameters allowing the charac¬ 
terization of planetary surface materials (e.g., grain size, morphology, internal 
structure and surface roughness). 

More recent developments are available. First, the version from |Hapke[ |2QQ2 
includes: (i) a more accurate analytic approximation for isotropic scatterers, 
(ii) a better estimation of the bidirectional reflectance when the scatterers are 
anisotropic, and (iii) the incorporation of coherent backscattering. Second, the 
version from |Hapke[ |2QQ8 treats the porosity 


Following previous studies (Johnson et ah, 2QQ6b|a Jehl et al. 2008 


Beck 


et al. 2012 Fernando et al. 2013), we decided to use the following standard 
expression of Hapke’s 1993 in order to be coherent with older analysis. In 
addition, more recent developments are not fully validated with experimental 
data. We remind here the main expression: 


to D \ ^ 

r{0o,0,g) = 


471 (lloe + lie) 

Using the following quantities: 


{[1 + B{g)] P{g) + Hiiioe)H{iie) - 1} 5(^0, e, g) 

( 1 ) 


• 6>o, 6>, and g: incidence, emergence and phase angles. The whole geometry 
quantities are noted O = O^g). Terms /ioe and are the cosine of the 

equivalent incidence and emergence angles, in the case of a rough surface, 
as defined in (Hapke 1993[ ). We note (f) as the azimuth angle. 
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• u; (0 < cj < 1): single scattering albedo. It represents the fraction of 
scattered to incident radiation by a single particle (also noted w). 


• P{g): particle scattering phase function. It characterizes the angular dis¬ 
tribution of energy for an average particle. We used the empirical 2-term 
Henyey-Greenstein function (Henyey and Greenstein 1941 1) (h ereafter re 
ferred to as HG2) for studying planetary surfaces (Gord et al.l|2003l |Jehl 


et al.[ 2008} [Johnson et al.[|2QQ6b|a||Souchon et al. 201 1| Beck et aH|2Q12| 


Pommerol et al.[|2Q13 ): 


P{g) = {i-c) 


+ c- 


(1 + 26cos (g) + (1 — 26cos (g) + 


( 2 ) 


The HG2 function depends on two parameters: the asymmetry parameter 
6 (0 < 6 < 1) which characterizes the anisotropy of the scattering lobe 
(from 6 = 0, which corresponds to the isotropic case, to 6 = 1, which 
corresponds to a particle which diffuses light in a single direction) and 
the backscattering fraction c (0 < c < 1) which characterizes the main 
direction of the diffusion (c < 0.5 corresponding to forward scattering and 
c > 0.5 corresponding to backward scattering). 


H{x): multiple scattering function. An analytical function for isotropic 


scatterers has been proposed (Hapke 1993) with a relative error to the 


exact values ( | Ghandrasekhar 1960) lower than 1%, leading to a relative 
error on the BRDF lower than 2% (Gheng and Domingue 2000). Defining 
^ = (1 — the multiple scattering function is: 


H{x) = \ l-[l-y]x 


^-y 

l + y 


1 -- 

2 


^-y 

l^y 


— X 


^-y 

l^y 


In 


1 X 


(3) 

A new expression dedicated to anisotropic scattering has been proposed 


(Hapke 2002). Nevertheless, Pommerol et al. (2013) have noticed that 


the use of the recent H expression leads to no significant changes over the 
previous expression. 

• B{g): opposition effect function. It describes the sharp increase of bright¬ 
ness around the zero phase angle often observed in the case of particulate 
media. Only the Shadow Hiding Opposition Effect (SHOE) is taken into 
account as follows (Hapke 1993): 


B{9) = 


Bo 


1 + ban(i) 


(4) 


The function depends on the parameters h and Bq (ranging from 0 to 1) 
which are the angular width and the amplitude of the opposition effect 
respectively. The Goherent Backscattering Opposition Effect (GBOE) is 
ignored in our case since the minimum phase angle is large {g > 10°). 
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S{0o, 0, g): macroscopic roughness factor. It describes the surface topogra¬ 
phy as a set of facets with a Gaussian distribution of tangent of the slopes, 
with mean slope angle noted 0 also called surface macroscopic roughness 
(Hapke 1993). This model of the surface roughness effect includes the 


partially shadowed area depending on the geometry and the bias on the 
effective incidence and emergence angles. The expression of S is given in 


Hapke (1993). 


2.2. Bayesian inversion 

Because the Hapke model is a non linear model and may not have a unique 
solution, one may use the Bayesian inversion framework based on the concept of 
the state of information which is characterized by a Probability Density Function 
(PDF) ( [Tarantola and Valette 1982). This approach has already been proposed 
for the Hapke model (Fernando et al. 2013). To infer the solution, this approach 


takes into account the initial state of information (a priori knowledge) on the 
parameters and the observations and applies the Bayes’ theorem to estimate the 
final state of information, called a posteriori. The key points of the concept and 
framework assumptions are presented in the following: 


Direct model and related quantities: 


— The direct model F estimates the simulated data d, from model pa¬ 
rameters m: 

d = F{m) (5) 

— The parameters m include uo^b^ d, Bq and h. The parameters m are 
in the parameter vector space M = [0,1]^[0,45], since all parameters 
are between 0 and 1, except 0 between 0 and 45° . 

— The simulated data d is the collection = r (O^) for all ith peculiar 
geometry noted as fti = (do, 0^g)i. The total number of geometries is 
noted Ng. The simulated data d are in the observation vector space 
since the BRDF is a positive quantity. 

— The model F{m) is the Hapke equation We consider that the 
geometries have very low uncertainties, which is the case for most 
data in planetary science cases. Thus, the parameters Vti are not 
estimated by the inversion. 

• Observation and other a priori information: 

— The actual observation is considered as prior information on data 
pD{d) in the observation space D. It is assumed to be a A^^-dimension 
gaussian PDF ^(o, C), with mean o and covariance matrix C. The 
values Oi are the observation for each geometry. The covariance ma¬ 
trix C is assumed here to be diagonal since each measurements at a 
given geometry is independent of the other geometric measurements. 
The diagonal elements Cu are ,..., , with being the stan¬ 

dard deviation. 


5 


















— The prior information on model parameters pM (^) in the parameter 
space M is independent to the data and corresponds to the state of 
null information if no information is available on the parameters. We 
consider an uniform PDF in their definition space M. Outside M, 
the PDF is null, avoiding unphysical solutions to appear. 

— The state of null information poid), representing the case when no 
information is available, is trivial in our case and represents the uni¬ 
form PDF in the parameters space M. 


• Solution of inverse problems and a posteriori information: 


— The posterior PDF in the model space as defined by Tarantola 


and Valette (1982) is: 


(Tm (m) = k Pm (m) L{m) (6) 

where fc is a constant and L{m) is the likelihood function 

£(m)= f (7) 

Jd PD\d) 

where 0{d\m) is the theoretical relationship of the PDF for d given 
m. We do not consider errors on the model itself, so 0{d \ m) = 
<5(F(m)). 

— The solution of the general inverse problem is given by the PDF 

o'Mi'm). The best way to represent crMi'm) is to plot the marginal 
PDF for one parameter j (see for instance fig. or the 

bivariate marginal PDF (TMij^j^mj') (see for instance fig.^. 

— The PDF can be described by statistic indicators such as mean values 
(expectation), standard deviations (covariance matrix), higher order 
statistics, etc. 


2.3. Monte Carlo Markov Chain (MCMC) to sample the inverse problem 

Since the Hapke model is non-linear, it is not possible to describe the pos¬ 
terior PDF cTMi^d) analytically. The solution provided by [Fernando et al. 


(2013) was to sample the final solution using a Monte Carlo approach using 


the Metropolis rule to built a Monte Carlo Markov chain (MCMC) (Mosegaard 


and Tarantola 1995). It consists in random testing of a candidate simulation 


over the parameter space and keeping only some solution in order to follow the 
a priori information. After a sufficient number of steps, the chain corresponds 
to the desired distribution, independent of the initialization. We used a very 
conservative number of 500 steps. 

The MCMC has been set to contain 500 sampling points. It may be not 
sufficient to have a smooth marginal PDF but more iterations are easy to com¬ 
pute using the same algorithm using more computation time. For the largest 
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shape in the bivariate PDF, we used 5000 sampling points instead. Let note the 
MCMC, sampling the final solution as: 

/ = 1,500 (8) 

The corresponding MCMC sampling of (TD{d) is estimated using: 


{fj, = ; = 1,500 


(9) 


et al.l 2003) and Particle Swarm Optimization (Beck et ah 2012 Pommerol 


Please note that this method is somehow similar to genetic algorit hms ([Cord 


et al.[ 2013), already used to solve this inverse problem empirically. Neverthe¬ 


less, those previous methods, although much faster, are only heuristic since no 
convergence proof of the algorithms have been proposed. This is not the case 
for the Metropolis rule, proposed here. 

The ability to rapidly find the “best solution” of heuristic method is very 
convincing but these methods are not able to estimate the uncertainties. In our 
case, the posterior PDF of a retrieved parameter is not necessarily a gaussian 
distribution which is the main advantage of the Bayesian method. Nevertheless, 
the Metropolis method can be very slow, especially when the solution is well 
constrained (i.e., a posteriori PDF with a very low standard deviation). 

In our case, all parameters have uniform physical prior distribution. If the 
solution has an uncertainty of 10% of the complete physical domain in M (for 
instance 0.5-0.6 for the parameter 5, which can vary from 0 to 1), the relevant 
subspace in the parameter space is only 0.1^ = 10“^ (for 6 parameters). It 
means that statistically, only 1 iteration over 10^ is kept in the MCMC and all 
other results are erased. To improve the computation time in this situation, we 


propose to use a fast Monte Carlo method described in section 2.5 


2.4‘ Description of the MCMC 

To describe the solution crM(^), in addition to simple histograms, additional 
statistical indicators on the MCMC are proposed, following [Fernando et al.| 

• The average value ifij (mathematical expectation) of each parameter j, 
and the estimated reflectance at each geometry i. 

• The covariance matrix Cm in the model space can be estimated from the 
MCMC. The dj standard deviation error bars on each parameter j are 
estimated from the covariance matrix elements Cmjj = cr|. 

• The non-uniformity criterion k. The parameters m are constrained if 
their marginal posterior PDF differs from the prior state of information 
(i.e., uniform distribution in our case). In order to distinguish if a given 
parameter is constrained we use the non-uniformity criterion k. 

Central moments iin (such as the variance /i 2 at order two) are commonly 
used for statistical purpose while cumulants kn have the advantage to 
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present unbiased statistical estimator for all orders ((Fisher 193Q| )). The 
first four cumulants /ci, /c 2 , /cs, of a uniform PDF are equal to 
0, — We proposed the parameter k to estimate the non-uniformity, 
defined as: 

h- I ^2- ^ ks k^- 
^ ’ 1 

60 120 


k = max 


1 

120 


12 


( 10 ) 


We perform a numerical test of 10,000 uniform random vectors of 500 
samples (identical to the MCMC) and compute k for each vectors. Since 
the maximum is k=0A7 for the most extreme vector, we propose to con¬ 
sider non-uniform PDF for k > 0.5, which is true with a probability higher 
than 99.99 %. For the inversion purpose, since the a priori PDF on the 
parameters are uniform, if the results of the inversion on one parameter 
has k < 0.5, we conclude that this parameters is not constrained by the 
observations. 


2.5. Fast MCMC 

The naive but accurate Monte Carlo Metropolis rule described in the previ¬ 
ous section may be not applicable in the case of long computation time due to, 
either a large number of sampled geometries, or a small data uncertainty. In the 
first case, the computation time is large due to the time required to compute 
one direct model (one candidate model). In the second case, the computation 
time is large because the algorithm tests numerous non relevant parameter set 
solutions (within the model space M but far away from the actual solution). 

In order to speed up the computation time, it is possible to use an adaptive 
Metropolis algorithm (Haario et al. 2001) but the solution (a posteriori PDF 
<JM{'m)) has to be close to a gaussian distribution. This algorithm uses the 
last Markov Chain to estimate a more relevant prior information pM{m) = 
Q{m^Cm)i recursively closer to the actual solution. If the solution is gaussian, 
the convergence of this method has been proven (Haario et al. 2QQ1| ). Also, 
this method is very convenient since it reduces very significantly the number 
of steps. In our tests, the speedup can reach a factor 100 without significant 
differences, when the final solution is well constrained. Table indicates an 
example of results on a full BRDF (see fig. 17 and 18). Differences between 
estimated parameters rh for MCMC and fastMCMC are always much lower than 
the estimated standard deviation a so we can consider as statistically equivalent. 


3. Synthetic tests 

We perform several synthetic observations in different conditions, to prop¬ 
agate the uncertainties from observations into the uncertainties on the Hapke 
parameters. 

The first goal is to determine favorable geometric conditions to accurately 
estimate the parameters for future spaceborne, in situ and laboratory investi¬ 
gations (section 1^. We will study two difficult cases : EPF observation and 
one single image / disk resolved image. The case of a full BRDF is not relevant 

























algorithm 

b 

c 

0 

UJ 

Bo 

h 

m 

MCMC 

0.45 

0.45 

7.42 

0.89 

0.50 

0.48 


fastMCMC 

0.46 

0.46 

6.13 

0.89 

0.47 

0.44 

a 

MCMC 

0.06 

0.08 

4.30 

0.02 

0.28 

0.29 


fastMCMC 

0.04 

0.06 

4.57 

0.02 

0.30 

0.27 

k 

MCMC 

1.00 

1.00 

0.99 

1.00 

0.25 

0.10 


fastMCMC 

1.00 

1.00 

1.01 

1.00 

0.25 

0.44 


Table 1: Comparison between results of the retrieved parameters ca, 6, c, Bq and h using 
classical Metropolis MCMC and fast MCMC, on average value m, standard deviation a and 
non uniformity criteria k. The data is a full BRDF measurement at 96 geometrical conditions 
and using the following parameter set: b = 0.5, c = 0.5, 6 = 1°, a;=0.9, h=0, Bo=0, and 10% 
uncertainties as described in fig. [^and[^ . 


because the uncertainties are small for all parameters cj, 6, c, 0 with coherent 
values, as shown in the example in table The opposition effect parameters 
Bq and h can only be constrained for small phase angle measurements and are 
out of the scope of this article. Those parameters are often studied separately 
in the literature. 

The second goal is to estimate if the Hapke hockey stick relationship could 
be due to a non-linear effect on the inversion (section]^. 

For each test, we compute a perfect model in REFF (REFlectance Eactor) 
unit at known geometry = {0o^0^g)i using eq. and known parameter mj. 
The reflectance in REEE unit is REFF = tt • r{0o^0^ g)/cos{0o). We model 
the uncertainties on the measurement as a gaussian function, independent from 
each geometry. The standard deviation level at geometry i is set to 10 % of 
the observed reflectance Oi in all the numerical tests, except when specified: 


10 


( 11 ) 


This value may be an upper limit for some spaceborne/laboratory instrumental 
uncertainties but it shall give an upper limit of the final uncertainties on the 
parameters. Also, taking all sources of error (including the atmosphere correc¬ 
tion), a noise level of 10% is realistic (Ceamanos et al.l |2013l lEernando et al. 


|2013[ ) . Eor the case of CRISM data, the reflectance error at each geometry were 


estimated at = r^/SO for instance (Eernando et al. 2013) 


4. Uncertainties propagation 

4.F Emission Phase Function (EPF) 

An Emission Phase Eunction (EPE) is a special configuration of observation 
with one particular incidence direction (incidence angle Oq) and a collection of 
emergence angle (emergence angles 0) along one single azimuthal plane. Table 
summarizes the different EPE conditions used. The EPE reflectance set, usually 
represented as a function of the phase angle g, is also called a photometric curve 
(Eigure[^ black curve). 
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Angular configurations 

6 »o (deg.) 

30 ; 40 ; 50 ; 60 ; 70 ; 80 

0 (deg.) [case # 1 ] 

[70, 50, 30, 10, -10, -30, -50, -70] 

0 (deg.) [case # 2 ] 

[70, 64, 58, 52, 47, 0, -47, -52, -58, -64, -70] 

0 (deg.) [case #3] 

[80, 70, 60, 50, 40, 30, 20, 10, 0, -10, -20, -30, -40, -50, -60, -70, -80] 

; ^2 } (deg.) 

{0 ; 180} ; {30 ;150} ; (60 ;120} ; (90 ; 90} 


Table 2: Angular configurations used in the EPF synthetic tests.We tested: (i) 6 incidence 
angles (h) 3 sets of emission angle 6 . The case #1 represents a standard EPF, case#2 
corresponds to a CRISM-like EPF observations and the case #3 corresponds to a favorable 
condition with a broad emission angle sampling typically used in laboratory investigations, 
(iii) 4 sets of azimuth angle Lp {^pi: azimuth angle in the illumination direction, (p 2 ’ azimuth 
angle in the opposite illumination direction) 


4 . 1 . 1 . Uncertainties in one EPF example 

We simulate a standard EPF observations (Table case 7 ^ 1 ) with a inci¬ 
dence angle Oq = 60° along the azimuthal plane {ipi ; (f 2 }={30; 210 } resulting 
in a phase angle range from 29° to 122° , and using the following model param¬ 
eter set: uj = 0.9, b = 0.8, c = 0.1, ^ = 15°, 5o = 0 and h = 0 corresponding 
to a bright material with a narrow forward scattering behavior and rough topo¬ 
graphic surface. This surface corresponds to granular soil with small grain size 
and round shape. A similar set of parameters have been observed in the labora¬ 
tory measurements of olivine at 700 nm (Souchon et al. 2011). Then, we invert 
the synthetic dataset with an uncertainty level set at 10 % of the reflectance (one 
standard deviation). We examine the final solution estimated from the last 500 
iterations of the Markov Chain. One have to remind that the discrepancies be¬ 
tween the solution and the initial data are not due to the retrieval method itself, 
but by the lack of information in the available data (poor geometric sampling, 
large uncertainties). 

Figurepresents the synthetic data and the fit of the 500 sampled solutions. 
The solutions are fitting the synthetic data with the expected tolerance (95% 
of the fits inside the 2 a data error). 

Figure shows the histogram of the Hapke parameters estimated from the 
500 sampled solutions. The opposition effect parameters BO and h present a 
flat histogram and have a non-uniformity criterion k < 0.5, suggesting that 
both parameters are not constrained. These results are directly related to the 
lack of observations at phase angle lower than 20 ° to observe the opposition 
effect. The single scattering albedo histogram shows a double peak around the 
expected value (c(;=0.9), showing the effect of the Hapke model non linearity. 

This example shows that the Bayesian inversion is able to identify several 
possible solutions due to data uncertainty and/or geometric sampling. The 
apparition of the double peaks on the parameter uj has been identified to chase 


the limitations of geometric diversity and/or reflectance uncertainty (Fernando 


et al. 2013). Since, usually uj is the best-constrained parameter in photometric 


modeling, the double peak in uj is also an indicator of low-constrain on other 
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Figure 1: Example of a result on a synthetic observation with 10% uncertainty. The black 
curve represents the synthetic data with one and two standard deviations. Light color curves 
represent 500 sampled solutions from the Monte Carlo Markov Chain. 


Hapke parameters. The parameters c and 0 show a broad PDF with a maximum 
close to the expected solution. Interestingly, the parameter b has a PDF that is 
not peaking to the expected solution, most probably because the phase range 
29° -122° at 10% uncertainties is not sufficient to constrain it. 

Figure [^presents the bivariate histogram for couples of parameters, permit¬ 
ting to study the combined effects of two parameters. It is clearly demonstrated 
that all parameters are correlated. For instance, the 6 vs c histogram clearly 
shows a “U” shape covering a large part of the model space M, but with a 
strong correlation, excluding medium b and strong c solutions. The single scat¬ 
tering albedo uj is better estimated with low c values, but higher c values are 
compensated by slightly lower cj, demonstrating the complex correlation of uj 
and c. 

This test shows that the Bayesian inversion is able to sample complex solu¬ 
tions in model parameters M which is not the case for classic inversion proce¬ 
dures based on minimization techniques. It also demonstrates that the maxi¬ 
mum likelihood may be a wrong indicator of the whole solution and that uncer¬ 
tainties may have a very complex shape due to correlation between parameters. 

4 . 1 . 2 . Favorable geometric sampling and uncertainties 

The objective of this sub-section is to evaluate the favorable geometric con¬ 
ditions in order to better estimate the Hapke photometric parameters by sim¬ 
ulating various surface scattering properties. These tests allow to evaluate the 
influence of the EPF geometric sampling in the uncertainties of the retrieved 
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Figure 2: Probability Density Function (PDF) of the solution, for each parameter of the 
Hapke model on the same example as fig. ^ Each plot represents the histogram of the 500 
acceptable solution from the Monte Carlo Markov Chain. The color line represent the initial 
parameters set. 


model parameters. These synthetic results should give constraints for current 
spaceborne data analyses (e.g., CRISM/MRO, HRSC/MRO, OMEGA/MEx for 
Mars, VIMS/Cassini for Titan, VIRTIS/ Rosetta for 67P/Churyumov-Gerasimenko, 
LROG/LRO for the Moon, etc) and future instrumental and laboratory inves¬ 
tigations. The following tests are done only for realistic soil cases, with hjc in 
the Hapke hockey-stick relationship. The opposition effect is not studied here 
because it is usually investigated separately from the global photometric data, 
this effect being only rarely observable in spaceborne data. 

We perform synthetic tests, using various surface scattering properties in 
order to cover the range of properties which can be observed in natural envi¬ 
ronments (i.e., different Hapke parameter sets): 

• 2 single scattering albedos (cj): 0.3; 0.9 

• 6 parameters (6 and c couples) taken along the hockey stick relation: 
Ll(6=0.3/c=1.0), L2(0.3/0.8), L3(0.3/0.5), L4(0.5/0.2), L5(0.8/0.1), L6(0.8/0.1) 

• 2 macroscopic roughnesses (^): 0° ; 15° 

• Opposition effect parameters are set to ignore this effect: Bo=0 and h=0. 

We perform the synthetic tests under various geometric configurations summa¬ 
rized in Table 6 incidence angles, 4 azimuthal modes and 2 cases of emission 
angles samplings where the case ^2 to a GRISM-like EPE observation with a 
poor-sampling and a narrow emission angle range (11 configurations) and the 
case 7^3 corresponds to a EPE with a well-sampling and broad emergence an¬ 
gle range (17 configurations). In what follows, we focus on the results of the 
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Figure 3: Probability Density Function (PDF) of the solution^or each couple of constrained 
parameters (cn, 6, c, 0) on the same example as fig. and The black/white diamonds 
represent the average of the PDF. The green circles represent the expected values for each 
parameter. 


13 






















case 7^2 as for example but all the results of the case 7^8 are presented in the 
Supplementary Materials. 

All 576 configurations are summarized in a single graph by computing the 
difference between initial parameter value of the synthetic test rrij and the esti¬ 
mated one from the Bayesian inversion fhj (with j the index of parameter and 
I the index of the Markov Chain): 


{/\mj}^ = {7flj}^-mj ( 12 ) 

This quantity represents the deviation between the initial “true” parameters 
and the state of information present in the EPF data, given the geometric 
sampling and uncertainties. The histogram of represents the PDF of 

uncertainties. By tracing marginal PDF, it is possible to chase the effect of each 
parameter on the uncertainties. One has to remind that the deviation between 
the solution and the initial value are not due to the retrieval method itself, but 
to the lack of information in the available dataset. 

Figures 1^ to present the deviation for one studied parameter each: geo¬ 
metric configurations: incidence angles (fig. [^, azimuthal modes (fig. [^, model 
parameter: c and b couple (LI to L 6 ) (fig. [^, uj (fig. and 0 (fig. by averag¬ 
ing the effects of all other parameters. The results show that the distributions 
of {Amj}i are always maximum on 0 , showing that the maximum likelihood 
is globally unbiased. It demonstrates that a single EPF is enough to estimate 
accurate model parameters in most of the cases. 

Nevertheless, the pick and the tails shape of the PDF of uncertainties are 
not identical for all marginal probabilities, showing that the CRISM-like EPF 
sampling can be optimized for some geometries, or for some surface parameters. 

One has to note that the roughness parameter is tested for the case ^ = 0 . 
Since 6 > > 0 by definition, the deviation cannot be symmetrical ({Am^}^ > 0 ), 
this significantly biases the deviation PDF. 

In detail, we can take the following conclusions: 

• The single scattering albedo parameter uj presents the narrower deviation 
PDF in comparison to other Hapke parameters (fig. [^to|^). Thus uu is 
the best constrained parameter on EPF measurement. 

• The influence of the geometric sampling on the deviation is very important: 
both incidence angle (fig. and azimuthal mode (fig. significantly 
change the global uncertainties: 

— Incidence has a strong effect on the estimation of surface roughness: 
the larger the incidence, the better the estimation of 0 . 

— Azimuth has a strong effect on 5, c and uj\ the closer to principal 
plane, the better the estimation. This effect is also present on 0 but 
more moderately. 

— We can interpret this behavior by the diversity of phase angles neces¬ 
sary to constraint the photometric parameters, both at low {g < 30°) 
and high {g > 100 °) ranges. Observations acquired at high incidence 
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angles and/or close to the principal plane allow to have the broadest 
phase angle ranges corresponding to the most favorable conditions. 
A departure from principal plane as low as 30° leads to significant 
increase of the uncertainties. 


• The influence of the model parameter values on deviation play a second 
order role: 


— The phase function parameters seem to have a strong influence on 
the retrieved parameters. In particular, the parameters 6, c and uj 
are significantly better estimated for surface materials with a narrow 
forward behavior (e.g., L6) (fig. |^. 

— The single scattering albedo uj has a moderate effect on the estimation 
of the parameters c and uj\ (i) the estimation of parameter uj is better 
when the albedo is high, (ii) the estimation of the parameter c is 
better when the albedo is lower (fig. [^. 

— The macroscopic roughness 0 has a small effect on the estimation 
of the parameters 6, c and uo: the parameter estimations are better 
when the parameter 0 is lower, (fig. [^. 


• Those conclusions are also valid for the case of favorable EPF (case 7^3, 
presented in the Supplementary Materials), except that the uncertainties 
are slightly reduced in this case. 


A similar study has been proposed on one experimental dataset using the genetic 
algorithm (Souchon et al. 2011). The authors recommended a “regular coverage 
of the bidirectional space in incidence, emission, azimuth, and consequently 
phase angles” and showed that “reliable photometric estimates can be produced 
with a limited set of angular configurations (on the order of a few tens)”. Our 
conclusions demonstrate that this experimental work can be extended even in 
the principal plane when large phase angle range are available. 


4.2. One single observation of a rough surfaee / disk resolved image 

One single image of a known rough surface has very limited angle geometry. 
Such observational condition corresponds to the case where each pixel of the 
image is assimilated to a facet with known orientation. We assume that the 
surface properties are spatially homogeneous over a large amount of pixel in 
order to estimate the surface properties by combining several adjacent pixels. 
This case is equivalent to a disk resolved measurement of a planetary body 
assumed to be homogeneous in surface properties, at one single phase image. 
We study this very difficult condition, in order to estimate if it is possible to 
constrain the photometric parameters. 

For each facet (or equivalently for each image pixel), the incident and emer¬ 
gent rays are parallel, with the same phase angle. However, since each facet 
has its own orientation, there is a variability of local incidence, emergence and 
azimuth angles. We propose to use a typical observation condition of incidence 
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PDF of the parameter b PDF of the parameter c 




Figure 4: Difference between true parameter and the Bayesian solution (see eq. |12| >, for 6 
incidence: 30° , 40° , 50° , 60° , 70° , 80° . All other parameters are included. 



Figure 5: Difference between true parameter and the Bayesian solution (see eq. |12| ), for 4 
azimuth: 0° , 30° , 60° , 90° . All other parameters are included. 
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PDF of the parameter b PDF of the parameter c 



PDF of the parameter theta-bar PDF of the parameter w 



Figure 6: Difference between true parameter and the Bayesian solution (see eq. |12| ), for 6 con¬ 
figurations of phase function parameters in the hockey stick: #1(6=0.3/c=1.0), #2(0.3/0.8), 
#3(0.3/0.5), #4(0.5/0.2), #5(0.8/0.1), #6(0.8/1.0) All other parameters are included. 


PDF of the parameter b PDF of the parameter c 



Figure 7: Difference between true parameter and the Bayesian solution (see eq. |12| >, for 2 
single scattering albedo uo\ 0.3, 0.9. All other parameters are included. 
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Figure 8: Difference between true parameter and the Bayesian solution (see eq. |12| ), for 2 
macroscopic roughness 6 \ 0° , 15° All other parameters are included. 


Oq = 40° , emergence 6> = 0° on 24 facets with 6 slopes 6>s = 5° , 10° , 20° , 30° 
, 40° , 50° and 4 azimuth of slope = 0°, 90° , 180° , 270° . Each facet is 
thus defined by a couple of slope and azimuth angles 0s). The phase angle 
is always g = 40°, but the local incidence Oq on the facet varies from 0° to 90° 
, the local emergence 0 from 5° to 50° and the local azimuth 0 from 0° to 180° 
. The relationship between slopes, azimuth of slopes and local incidence and 
emergence are the following Hapke (1993): 


cos ^0 = cos ^0 cos Os + sin Oq sin Og sin 0^ 


(13) 


COS 0 = cos 0 cos Os + sin 0 sin Og sin 0^ (14) 

We test uj = 0.9, 0 = 1^ and three values of the h/c in the hockey stick : 
#l(6=0.3/c=1.0), #4(0.5/0.2), #6(1.0/0.1). 

The results are plotted in figure in terms of bi-variate histograms describ¬ 
ing b vs c, and uj vs 0. The parameter uj seems to be well constrained in all cases 
but the parameter 0 is only constrained in the case of a strong narrow forward 
scattering. The particulate phase function parameters b is only constrained in 
the extreme cases of the hockey stick (cases #1 and #6). The parameter c is 
not constrained. 

These results indicate that information of the surface properties can be re¬ 
trieved even in the case of one single observation with known geometry for each 
facet. Especially uo can be retrieved with small uncertainties (la uncertainties 
< 0.1 ), but also 6, c and in some extent 0 , in the case of very extreme phase 
function. 
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Figure 9: Probability Density Function of the solution, for each couple of constrained pa¬ 
rameters (ca, 6, c, 6 ) derived from the simulation of a single observation of a known rough 
surface with 24 facets (see the text for more detail). At the left case #l(6=0.3/c=1.0). In the 
middle, case #4(0.5/0.2). At right case #6(1.0/0.1). The black/white diamonds represent 
the average calculated from of the PDF. The green circles represent the expected values for 
each parameter. 
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5. Possible origin of the hockey stick relationship 


We exhibit here particular conditions of Emission Phase Function (EPF) 
and Bidirectional Reflectance Distribution Function (BRDF) that could be a 
possible bias at the origin of the hockey stick relationship, so we focus in this 
section on the parameters b and c. We set cj=0.9, 0 = 1° , h=0, Bo=0 and data 
uncertainties at a level of 10% as previously defined. T he res ults are expected 
to be only weakly dependent of uj and 0 (see section 4.1.2). These surface 
properties correspond to granular soil with small grain size. Similar parameters 
have been observed in sulfate terrain on Mars (Johnson et al. 2006a), and in 
various samples (Souchon et al. 2011 [ ), such basalt including feldspar grains. 


pyroclastic grains and olivine grains. 

We test here all configurations of 6/c covering the whoie parameter spaces 
from 0 to 1 and see if the resuiting uncertainties can bias the interpretation or 
not. We insist on the fact that materiais with both h> 0.5 and c > 0.5 have 
never been observed. Nevertheiess, we tested these configurations in order to 
study the potentiai bias when anaiyzing such case. 


5.1. Principal plane EPF 

Figure [^presents the resuits for a standard principai piane EPF observation 
(Tabie 1^ case 7 ^ 1 ) (i.e. poor-sampling of emission angles, mostly the case for 
spaceborne instruments like CRISM). It shows that low b (< 0.5) imply poor 
constraint on c and very large uncertainties on b. For 6 > 0.5 and c > 0.5, b has 
small uncertainties and c has medium uncertainties. 

Figure prlpresents the results for a favorable principal plane EPF observation 
(i.e., well sampling of emission angles which can be obtained at laboratory, on in 
situ in very favorable conditions). The solution is clearly better constrained than 
for the standard principal plane (see fig. [^. It demonstrates that even at 10% 
data uncertainties, the correct behavior (backward/forward and narrow/broad) 
can be retrieved from a single EPF observation if the observation is taken in or 
close to the principal plane. 


5.2. Effect of azimuth on EPF 

Figures pT| p!3l and p!5] present the results of the favorable EPF with an 
azimuthal plane, respectively {0°,180°}, {30°, 210°}, {45°, 225°}, {60°, 240°}, 
{90°, 270°}. As expected, the information included in one EPF observation 
to constrain the parameters b and c is decreasing by increasing the azimuthal 
plane angle. At an azimuthal plane angle of 0° (principal plane, shown in fig. 
pr] ) , the solutions are well constrained when the parameter b value is greater or 
equal 0.5 as observed in section 5.1 At an azimuthal plane angle as low as 30° 


(Figure [T^, the solution is not well constrained in most cases, except for the 
cases when 6=0.5 and the case when 6=0.9 and 6=0.1. At an azimuthal plane 
angle of 90° (Figure [^, only the parameter 6 can be qualitatively estimated 
(i.e., by discriminating the broad and the narrow scattering), the parameter c is 
unconstrained (i.e., all the solutions for the parameter c are possible), because 
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Figure 10: Probability density function of the b and c parameters using a standard principal 
plane EPF observation (8 emergence configurations: 70° , 50° , 30° , 10° , -10° , -30° 
, -50° , -70° ), incidence=40° , azimuth={0,180°} and using the following parameter set: 
00=0.9, ^ = 1° , h=0, Bo=0, data uncertainties 10% and 9 combinations of 5=0.1/0.5/0.9 
and c=0.1/0.5/0.9. The black/white diamonds represent the average of the PDF. The green 
circles represent the expected values of the parameters b and c. 
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Figure 11: Probability density function of the b and c parameters using a favorable principal 
plane EPF observation (17 emergence configurations: 80° , 70° , 60° , 50° , 40° , 30° , 20° , 10° 
, 0° , -10° , -20° , -30° , -40° , -50° , -60° , -70°_, -80° ), incidence=75° , azimuth={0,180°} 
and using the following parameter set: a;=0.9, ^=1° , h=0, Bo=0, uncertainties 10% and 9 
combinations of 6=0.1/0.5/0.9 and c=0.1/0.5/0.9. The black/white diamonds represent the 
average of the PDF. The green circles represent the expected values of the parameters 6 and 
c. 
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no phase angles greater than 90° are available in the data, very important to 
distinguish the forward and the backward scattering direction. 

Interestingly, the behavior of a high b value coupled with moderate to high c 
value differs from the other cases: a “U” shape solution is expressed only in this 
quadrangle. If the tool used for the inversion is based on a root mean square 
minimization, it clearly depends on the initialization. If the initialization is in 
the hockey stick shape, it would converge to a fake local maximum. This effect 
may be a significant bias on the estimation of the parameters b and c, leading 
to an artificial hockey stick effect, present for EPF data with azimuthal plane 
larger than 30° and inappropriate inversion method. 


5.3. Effect of noise level on EPF 

In order to test the noise level on the retrieved parameters, we vary the noise 
level (eq. pT| , from 50% to 1%. First, figureclearly shows that with 50% data 
uncertainty, all solutions in the parameters b and c spaces are possible, thus the 
parameters are not constrained due to the large data uncertainty. With 10% 
data uncertainty, the solution is only restricted to the “U” shape leading to the 
artificial hockey-stick like shape. The b and c couple can be correctly evaluated 
with a data uncertainty lower than 5%. Such uncertainty level can often be 


obtained in laboratory measurements (Pommerol et al. 2013 Johnson et al 


2013) but sometimes also in spaceborne measurements (Fernando et al. 2013) 


Finally with 1% data uncertainty, the solution is well constrained close to the 
expected solutions. 

The artificial hockey stick effect discussed in the previous section is only 
present for uncertainties larger than 5-10%. 


5.4. Full BRDF 

We test the BRDF observation, sampled 48 times at 2 incidence angles 40° 
and 60° , 8 emergence angles: 70° , 50° , 30° , 10° , -10° , -30° , -50° , -70° 
and 3 azimuth angles: 0° , 45° and 90° . This sampling corresponds to the 
combination of standard EPF (table case 7 ^ 2 ) for 3 azimuthal plane angles, 
0° , 45° , 90° and two incidence angles, allowing high diversity of geometries. It 
represents typical laboratory measurements or the best expected combination 
of single EPF on the same site. Such combination of EPF have already been 


et al. 


proposed in the literature (Jehl et al. 2008 Fernando et al. 2013 Fernando 


2015). 


The results, presented in fig. show that the BRDF configuration contains 
enough information to constrain the parameters b and c, in comparison with 
single EPF (fig. [Tq| ) with lower uncertainties. Interestingly, the favorable EPF 
condition (fig. |11| ) is able to better constrain the phase function parameters 
than the standard BRDF presented here. This effect is most probably due to 
the maximum phase angle that is limited to 130° for the BRDF considered here 
but goes to 155° for the favorable EPF. 

A full BRDF observation at 10% uncertainties is able to constrain the pa¬ 
rameters b and c, and should reduce the contribution to the hockey stick artifact. 
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Figure 12: Same as fig. pT]but with azimuth={30°, 210°}. 
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Figure 13: Same as fig. pT]but with azimuth={45°, 225°}. 
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Figure 14: Same as fig. pT]but with azimuth={60°, 240°}. 
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Figure 15: Same as fig. |ll| but with azimuth={90°, 270°}. 
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Figure 16: Probability density function of the parameter b and c using a favorable EPF 
observation (17 emergence configurations: 80° , 70° , 60° , 50° , 40° , 30° , 20° , 10° , 0° , 
-10° , -20° , -30° , -40° , -50° , -60° , -70° , -80° ), incidence=75° , azimuth={45°, 225°} and 
using the following parameter set: u;=0.9, ^ = 1° , h=0, Bo=0, 5=0.9 and c=0.6. Different 
uncertainties values are tested: 50%, 10%, 5%, 1%. The black/white diamonds represent the 
average of the PDF. 


6. Discussions and Conclusion 

We proposed a rigorous inversion scheme to estimate Hapke model param¬ 
eter from measurements, using Bayesian Monte Carlo strategy. The typical 
computation time on a 2.5 GHz Intel Core 15, 8Go RAM, is one minute for a 
single EPF but can reach one hour for a BRDF. In order to speed up the con¬ 
vergence, we developed a Fast Monte Carlo strategy, reducing the computation 
time by a factor of 100. This strategy is only suitable on full BRDF or favorable 
EPF, when gaussian like solutions are expected. 

We explored various conditions on synthetic examples: EPF type observa¬ 
tions, one single observation, BRDF type observations in order to study the 
propagation of errors from measurements to the parameter space with 10% un¬ 
certainties as a nominal condition. The major conclusions of this work are: 

• Non-linearities in the Hapke model are important for EPF type measure¬ 
ments leading to potential multiple solutions, at least with data uncer¬ 
tainties larger than 5% and large azimuthal plane angle (> 30°). 

• Azimuthal plane in a EPF observation is the most important parameter to 
constrain the photometric parameters: the closer to the principal plane, 
the best the results. A departure of only 30° in azimuthal plane leads to 
significant increase of uncertainty. 

• Incidence angle is very important to constrain the parameters in a EPF, 
especially the roughness parameter 0 . A recommendation for laboratory or 
spaceborne observation is to sample the highest incidence angle possible. 

• One single EPF type observation with very favorable conditions (i.e., prin¬ 
cipal plane, incidence at 75° , emergence angle sampling up to 80° ) is 
enough to constrain cj, 6, c, and 0 parameters, even with a data uncer¬ 
tainty level of 10%. 


28 












Pdf (%) 


Pdf (%) 


Pdf (%) 



Figure 17: Same as fig. |10| but for a full BRDF observation including 48 geometries (2 
incidence angles: 40° and 60° , 8 emergence angles: 70° , 50° , 30° , 10° , -10° , -30° , 
-50° , -70° along 3 azimuth angles: 0° , 45° , 90° and using the following parameter set: 
a;=0.9, ^=1° , /i=0, 5o=0, data uncertainties 10% and 9 combinations of 6=0.1/0.5/0.9 and 
c=0.1/0.5/0.9. 
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Figure 18: 


Same as fig. ^|but using the Fast Monte Carlo method explained in section 
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• For data uncertainty less than 5%, the parameters can be estimated using 
single EPF under certain geometric configurations: close to the principal 
plane (azimuthal angle less than 45° ) and high incidence angles (greater 
than 50° ) leading to a broad phase angle range, containing low and high 
phase angles, to sufficiently describe the shape of the photometric curve. 

• In the case of one single observation, with each pixel considered as a facet 
with known geometry, uj can be retrieved with small uncertainties (la 
uncertainties < 0.1 ), but also 6, c and in some extent in the case of very 
extreme phase function. This case is equivalent to one disk resolved image 
of a planetary body, assumed to be homogeneous in surface properties. 

• Full BRDF observations allowing a high diversity of geometric sampling 
are the best configurations to constrain all the parameter set: cj, 6, c and 
0 . This geometric conditions can be easily reproduced in laboratory and 
can be obtained by combining different EPF observations at varied illu¬ 
mination conditions and/or varied azimuthal planes. Nevertheless, even 
BRDF measurements are limited by the phase range. Our results indicate 
that a favorable EPF with higher phase range than a BRDF is better to 
constraint the parameters. 

• A favorable EPF measurement (i.e., principal plane, incidence at 75° , 
emergence angle sampling up to 80° ) is better to constrain the parameters 
than a standard EPF (incidence at 40° and 60° , emergence up to 70° ), 
most probably due to better high phase angle sampling. 


• The hockey stick relationship on the 5 vs c diagram may be the result 
of the non-linearity of the Hapke model if the data are from a EPF type 
observation and the inversion strategy is based on simple mean square 
minimization. However, the Full BRDF type observations are not biased 
by the non-linearity. Because the data used in the Hapke (2012) synthesis 
are generally BRDF type observations (Hapke 2012), it is unlikely that 
the hockey stick relationship is an artifact from the inversion method. 
This confirms that surface material with strongly backward scattering 
with narrow lobe may not exist in the nature. 


Future work should include real laboratory spectra and datasets on planetary 
bodies with prioritization using the conclusion of this study. Also the wavelength 
dependance of all parameters should be addressed. Finally, the latest develop¬ 
ments of the Hapke model should be included within this inversion strategy in 
order to compare the actual properties of the granular material and retrieved 
photometric parameters, given precise uncertainties. 
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